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Diffraction and near-zero transmission of flexural phonons at graphene grain boundaries 
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Graphene grain boundaries are known to affect phonon transport and thermal conductivity, sug¬ 
gesting that they may be used to engineer the phononic properties of graphene. Here, the effect of 
two buckled grain boundaries on long-wavelength flexural acoustic phonons has been investigated 
as a function of angle of incidence using molecular dynamics. The flexural acoustic mode has been 
chosen due to its importance to thermal transport. It is found that the transmission through the 
boundaries is strongly suppressed for incidence angles close to 35°. Also, the grain boundaries are 
found to act as diffraction gratings for the phonons. 


I. INTRODUCTION 


Grain boundaries in graphene have been found to af¬ 
fect the mechanical, electronic and thermal properties of 
the materialiti^. The grain boundaries commonly consist of 
dislocations, in the form of pentagon-heptagon defect pairs, 
and cause out-of-plane buckling of the graphene sheet^ri^. 
Recent experimental studies show that dislocations can be 
introduced into pristine graphene using a focused electron 
beami^ii^di^— , suggesting the possibility of adjusting the 
properties of the material. 

The possibility of manipulating the properties of 
graphene could be particularly important in applications 
related to phononics and heat management22i2i, where con¬ 
trol of the vibrational properties and thermal conductivity 
of graphene is essential. The effect of grain boundaries 
on the thermal conductivity of graphene has previously 
been studied using both molecular dynamics and Greens 
function methods^^r.^^. However, out of these studies only 
Liu et ali^, who consider transport along the boundary, 
mention the influence of out-of-plane buckling. Also, these 
studies give no detailed insight into the scattering processes 
of specific phonon modes. 

In the present study, we investigate the scattering of long- 
wavelength flexural acoustic phonons at grain boundaries 
in graphene for several incidence angles using molecular 
dynamics (MD). This particluar phonon mode was chosen 
since it is believed to contribute significantly to the thermal 
conductivit y^^’^° . Two grain boundaries are considered in 
this paper, one with a misorientation angle of 9.4° and one 
with a misorientation angle of 17.9°. Both grain bound¬ 
aries display substantial out-of-plane buckling, with a pe¬ 
riodic variation in height along the grain boundary due to 
the distribution of defects. The boundaries are found to act 
as diffraction gratings for the phonons, and strongly sup¬ 
pressed transmission is also observed for specific angles. In 
particular, the transmission is as low as 4 % for incidence 
angles near 35° at both boundaries. 

A previous investigation limited to phonons normally in¬ 
cident on the grain boundary showed that the scattering 
was due almost entirely to the out-of-plane buckling of the 
boundary^. Based on this result a continuum mechani¬ 
cal model was constructed, where the grain boundary was 
modeled as a static out-of-plane displacement. The model 
showed good agreement with the MD results. Here, we ex¬ 
tend this continuum mechanical model to the case of non¬ 



Figure 1. (Color online) Symmetric tilt grain boundary with 
misorientation angle 9.4°, seen from the y direction (top) and a 
directon (bottom). Figure made using VMD^. 

normal angle of incidence in order to gain a qualitative 
understanding of the scattering mechanism. 


II. METHOD 

All MD simulations have been performed using the 
program package LAMMPS (large-scale atomic/molecular 
massively parallel simulator}^. The interaction be¬ 
tween carbon atoms has been modeled using the Tersoff 
potential^i^ with the potential parameters given by Lind¬ 
say and Broido^. This set of parameters has been chosen 
due to its improved description of acoustic phonon modes 
in graphene. The considered grain boundaries are symmet¬ 
ric tilt grain boundaries and consists of periodic arrays of 
pentagon-heptagon defects. The 9.4° grain boundary has 
a period of 1.5 nm in the y direction, parallel to the grain 
boundary (see Figure [1]), while the 17.9° boundary has a 
period of 2.4 nm. The grain boundaries have been con¬ 
structed using the method described in Ref. |3l|. For the 
9.4° boundary this results in a grain boundary buckling 
0.6 nm high and 1.7 nm wide. Due to the defect distribu¬ 
tion the buckling height varies periodically along the grain 
boundary with an amplitude of 0.06 nm. The 17.9° bound¬ 
ary has a buckling height of 1.5 nm and a buckling width 
of 5 nm, wih a variation of 0.1 nm along the boundary. 

To construct the phonon wavepackets we use the method 
of Kimmer et al^l. The displacement Uj of atom j is then 













Figure 2. (Color online) Symmetric tilt grain boundary with misorientation angle 17.9°, seen from the y direction (top) and 2 
directon (bottom). Figure made using VMD^. 


determined by 


Uj = Re 




( 1 ) 


where k = kxX + kyjj is a wavevector, e,k is a polarization 
vector for the considered phonon branch, Rj is the posi¬ 
tion vector of atom j and oj is the phonon frequency. The 
amplitudes Ok are calculated according to 

ttk = (2) 


In our previous study of phonon scattering at graphene 
grain boundaries a simple continuum mechanical model of 
the system was constructed in order to further confirm the 
results and to facilitate future studies of systems too large 
to model using MD^. The model built on the observa¬ 
tion that the main cause for scattering of long-wavelength 
phonons at the grain boundary is the buckling. Here, we 
have extended the previously used model from one to two 
dimensions for the case of the 9.4° boundary, and incorpo¬ 
rated the periodic height variation of the buckling. 

The equations of motion for the displacements are: 


where A is an amplitude and rj is the width of the 
wavepacket in the x direction (perpendicular to the grain 
boundary). The resulting wavepacket is localized in the 
X direction, centered around Rq in real space and around 
a wavevector ko = kxox + kyij in reciprocal space. All 
wavevectors k are required to be reciprocal lattice vectors 
of the simulation supercell. With periodic boundary condi¬ 
tions applied in the y direction (parallel to the grain bound¬ 
ary) this gives ky = 27r?7i/i®°, where m is an integer and 
is the size of the supercell in the y direction. 

The polarization vectors and dispersion relation a;(k) 
have been obtained from the dynamical matrix of the per¬ 
fect lattice using the General Utility Lattice Program, 
GULP— 1 ^, and the constants A and y have been set to 
0.013 and 5 nm, respectively. Since the focus of this study 
is long-wavelength phonons, the upper limit for |ko| = 

fc^o + ky has been set to 7 nm“^, which limits the possi¬ 
ble values of m and kxo- To extend the range of possible m 
values the size of the simulation supercells in the y direc¬ 
tion has been increased. For the 9.4° boundary it has been 
tripled, so that = 4.5 nm, while for the 17.9° boundary 
it has been doubled, giving = 4.8 nm. The supercells 
of the 9.4° and 17.9° boundaries are 260 and 400 nm long 
in the x direction, respectively. Fixed boundary conditions 
are applied in this direction and all atoms less than 10 nm 
from the supercell edge are held immobile. 


pil dxO'xx dydxy — 0 (3) 

pv dxC^xy ^y^yy — ^ (^) 

pw + - dx[(TxxdxW + axydyw] (5) 

^y\^xyk^x'^ 4 “ ^yy^y^] ~ 

where u is the displacement in a;, u is the displacement 
in y, ic is the out of plane displacement, p is the density, 
K is the bending rigidity and Oxx, <Xxy and Uyy are the 
components of the stress tensor. As in the previous study 
the grain boundary buckling has been included in the form 
of a static out-of-plane displacement. 

Finite-difference time-domain methods have been used to 
propagate wavepackets similar to the ones used in MD and 
to study scattering against the buckling. Results of these 
calculations can be directly compared to the MD simulation 
results. The details of the continuum mechanical model can 
be found in the Appendix. 


III. RESULTS 

The time evolution of the kinetic energy in both grains 
for a wavepacket with kxo = 4 nm and m = 2 interacting 
with the 9.4° boundary can be seen in Figure|31 Here, grain 
1 is defined as the grain in which the pulse is introduced. 
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Figure 3. (Color online) The fraction of the total kinetic energy 
in grain 1 (top) and grain 2 (bottom) as a function of time for 
a wavepacket with = 4 nm”*^ and m = 2 scattering at the 
9.4° boundary. 


and grain 2 is the other grain. Changes in the kinetic energy 
of the grains can be seen at two points. After 20 ps, the 
kinetic energy in grain 1 decreases to 73 % of the total 
kinetic energy while the kinetic energy of grain 2 increases 
to 27 %, indicating that the pulse has reached the grain 
boundary. The second change occurs at 60 ps, where the 
energy of grain 1 decreases further in two steps, first to 
60 % and then to 44 %. Between these two points the pulse 
has been reflected against the fixed boundary conditions, 
so that the steps at 60 ps mark the return of the scattered 
pulses to the grain boundary. 

The most surprising feature of Figure [3] is the stepwise 
change in energy beginning at 60 ps, which seems to indi¬ 
cate that there are two pulses arriving at the grain bound¬ 
ary about 5 ps apart. A closer examination of the scattered 
pulses shows that this is indeed the case. Figure 0] shows 
the intensity of the scattered pulses, normalized by the to¬ 
tal intensity, as a function of wavevectors and ky for 
t = 40 ps. Four peaks are seen, two with negative k^, cor¬ 
responding to reflected pulses, and two transmitted pulses 
with positive kx- The reflected pulses are labeled R1 and 
R2. R1 has kx = —4 nm“^ and ky = 2.8 nm“^, while R2 
occurs at kx = —4.7 nm“^ and ky = —1.4 nm“^. Similarly, 
the transmitted pulses T1 and T2 have fca, = 4, ky = 2.8 
nm“^ and kx = 4.7, ky = —1.4 nm“^, respectively. T1 
has the same wavevector as the incident pulse. Since the 
propagation velocity of the pulse depends on the value of 
kx , these two pulses will propagate with different velocities 
and thus give rise to the stepwise change in kinetic energy 
seen in Figure [S] 

The same phenomenon is observed at the 17.9° grain 
boundary. Figure [S] shows the normalized intensity after 
scattering for a pulse with kxQ = 4 nm“^ and m = 3. Four 
reflected peaks and three transmitted peaks can be seen. 


Figure 4. (Color online) Normalized intensity Ijhot after scat¬ 
tering at the 9.4° boundary (t = 40 ps) as function of kx and 
ky for a wavepacket with kxo = 4 nm“^ and m = 2. R1 and R2 
denote the reflected pulses, while T1 and T2 are the transmit¬ 
ted pulses. The dotted lines represent the values of ky allowed 
by the boundary conditions and the dashed circle indicates the 
points with ko = y'fcl + fcf equal to that of the incident pulse. 


For the reflected peaks, R1 occurs at kx = —4.0, ky = 
3.9, R2 at kx = —5.5, ky = 1.3, R3 at kx = —5.5, ky = 
— 1.3 and R4 at kx = —4.0, ky = —3.9 nm“^, while the 
transmitted peaks occur at kx = 4.0, ky = 3.9 (Tl), kx = 
5.5, ky = 1.3 (T2), and kx = 4.0, ky = -3.9 nm'^ (T3). 

Examination of the scattered pulses at both grain bound¬ 
aries reveal that the difference between the ky value for the 
incident pulse, k™, and the ky value for the scattered pulses, 
ky^, can be expressed as 


hsc _ T in _ 

S r ’ 

Ijy 


( 6 ) 


where n is an integer and Ly is the grain boundary pe¬ 
riod. The kx value of the scattered pulses, is given by 
momentum conservation: 


kx = ^{k^f - {kl^y- ( 7 ) 


This shows that the buckled, periodic grain boundaries act 
as diffraction gratings for long-wavelength flexural acoustic 
phonons. Grain boundaries functioning as diffraction grat¬ 
ings for phonons has previoulsy been used to model the 
behaviour of the thermal conductivity in ionic materials^. 

Unlike the previously described case where m = 0^, 
scattering into in-plane vibrational modes is negligible for 
all cases with m = 1 and m = 2. Some movement in the y 
direction is seen at the 9.4° boundary for m = 3. 

The transmission coefficient T is defined as 


T = 


^^grai„2) 
^tot ■ 


( 8 ) 
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Figure 5. (Color online) Normalized intensity I/hot after scat¬ 
tering at the 17.9° grain boundary as function of kx and ky for 
a wavepacket with k^o = 4 nm“^ and m = 3. R1 to R4 denote 
the reflected pulses, while Tl, T2 and T3 are the transmitted 
pulses. The dotted lines represent the values of ky allowed by the 
boundary conditions and the dashed circle indicates the points 
with ko = t^k% + ky equal to that of the incident pulse. 

where is the kinetic energy in grain 2, £1^°* is the 

total kinetic energy and the brackets represent a time av¬ 
erage over times between the first scattering at the grain 
boundary and the time when the first wavepackets reaches 
the edge of the supercell. Values of T for the 9.4° grain 
boundary for several values of at m = 1, 2 and 3 can 
be seen in Figure [6l For all values of m, the transmission 
increases with increasing k^o- The increase is monotonic for 
m = 3, while for m = 2 there is a small dip around kxo = 4 
nm“^ and for m = 1 there is a pronounced trough around 
kxo = 2 nm“^. Remarkably, the transmission for m = 1 
and kxo = 2 nm“^ nearly reaches zero, so that no part of 
the incident pulse is transmitted through the boundary. It 
can be noted that the dip in the curve for m = 2 and the 
trough for m = 1 occur at the same angle, but for different 
values of fcp. 

Figure [7] shows the dependence of T on kxo with m = 1, 2 
and 3 for the 17.9° boundary. As for the 9.4° boundary, the 
transmission increases with increasing kxo- Extremely low 
transmission is also observed at m = 1 and kxo = 2 nm“^, 
corresponding to an incidence angle of 33°. It is not clear 
whether there is a minimum at the same angle of incidence 
for m = 2, as in the 9.4° case, as the transmission is quite 
low also at slightly larger incidence angles. 

Figure [ 5 ] also contains transmission coefficients Tc ob¬ 
tained from the continuum mechanical model. The quali¬ 
tative agreement between the continuum mechanical model 
and the MD results is very good, as the continuum mechan¬ 
ical model clearly reproduces the general trend in the MD 
data of increasing transmission with increasing kxo- The 
two models agree particularly well for m = 3, although 
the continuum mechanical model overestimates the trans¬ 



Figure 6. (Color online) Transmission at the 9.4° boundary as 
a function of k^o for m = 1 (top), m = 2 (middle) and m = 3 
(bottom). The angle of incidence is indicated beside each data 
point. Open symbols and dashed lines represent results from 
the continuum mechanical model. 


mission at kxo = 5 nm“^. For m = 2, the dip around 
kxO = 4 nm“^ is reproduced but is wider than in the MD 
data, extending to kxo = 3 nm“^. The continuum mechan¬ 
ical model also overestimates the transmission at kxo = 5 
nm“^. Finally, for m = 1 the transmission obtained with 
the continuum mechanical model is higher than that ob¬ 
tained with MD over almost the entire interval. It also 
does not reproduce the trough at kxo = 2 nm“^, but does 
reach near-zero values for kxo = 1 nm“^. 

In addition to the transmission coefficient, the continuum 
mechanical model should reproduce the diffraction seen in 
MD. Figure [ 5 ] shows the intensity obtained from the con¬ 
tinuum mechanical model after scattering as function of kx 
and ky for kxo = 4 and m = 2, corresponding to the MD re¬ 
sults presented in Figure 0 ] It is clear that the same peaks 
appear, showing that diffraction occurs also in the contin¬ 
uum mechanical model. Compared to the MD results Tl 
appears to be underestimated and T2 overestimated, possi¬ 
bly due to that the model of the boundary buckling used in 
the continuum mechanical model does not reproduce the 
actual curvature of the grain boundary buckling in suffi¬ 
cient detail. 
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IV. CONCLUSION 



Figure 7. (Color online) Transmission at the 17.9° boundary as 
a function of k^o for m = 1 (top), m = 2 (middle) and m = 3 
(bottom). The angle of incidence is indicated beside each data 
point. 



Figure 8. (Color online) Normalized intensity after scattering as 
function of and ky for a wavepacket with k^o = 4 nm“^ and 
m — 2, from the continuum mechanical model. The dotted lines 
represent the values of ky allowed by the boundary conditions 
and the dashed circle indicates the points with ko = 
equal to that of the incident pulse. 


In summary, the effects of the angle of incidence on 
the scattering of long-wavelength flexural phonons against 
grain boundaries in graphene have been studied using 
molecular dynamics. The considered grain boundaries, two 
buckled symmetric tilt grain boundaries with misorienta- 
tion angles 9.4° and 17.9°, have been found to act as diffrac¬ 
tion gratings for long-wavelength flexural phonons. In ad¬ 
dition, near-zero transmission has been observed for angles 
near 35° and small wavevector magnitudes. A continuum 
mechanical model of the system containing the 9.4° bound¬ 
ary has been constructed and shown to qualitatively agree 
with the MD results, giving insights into the scattering 
mechanism and providing a starting point for studies of 
systems too large to be modeled atomistically. The pre¬ 
sented results improve our understanding of how phonons 
interact with grain boundaries in graphene and suggest that 
such defects could indeed be useful in manipulating the vi¬ 
brational properties of the material. 


Appendix: Continuum mechanical modeling 


In the continuum mechanical model the graphene sheet 
is described as a thin plate. The equations of motion for 
the displacements are 

pil dxCfxX dy(7xy — 0 , 
pv dxCTxy ^y^yy — 

pii) + kA'^w - dx[(JxxdxW + axydyw] 

^y\^xy^x^ 4 “ ^yy^y^] ~ 

where u is the displacement in x (perpendicular to the 
boundary), v is the displacement in y (parallel to the 
boundary), w is the out-of-plane displacement, k is the 
bending rigidity, p is the density and Uxx, cFxy and Uyy 
are the elements of the two-dimensional stress tensor. To 
model the grain boundary buckling a static out-of-plane 
displacement wo{x,y) is introduced. The introduction of 
this out-of-plane displacement gives rise to static displace¬ 
ments in the in-plane directions, so that the total displace¬ 
ments must be written 


(A.1) 

(A.2) 

(A.3) 


uix,y,t) = uo{x,y) + ui{x,y,t) (A.4) 

v{x,y,t) =vo{x,y) +vi{x,y,t) (A.5) 

w{x,y,t) =wo{x,y)+wi{x,y,t) (A.6) 


where ui(x,y,t), vi(x,y,t) and wi{x,y,t) are the time- 
dependent displacements. However, the displacements de¬ 
termine the stress tensor components through the relations 


^XX -4” 


^yy 


dxU ■ 


dxU 


(dxw) 


(dxw) 


+ X 


4- (A 4- 2p) 


dyV 


dyV 


{dyWY 


(A.7) 


{dyWf 


^xy - P 4- Oytl dxWdyW\ , 
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where A and fi are Lame parameters. Thus, the stress ten¬ 
sor elements can also be divided into a time-dependent term 
'^ij (b J = y) a time-independent term The two- 
dimensional stress tensor components , <7°^, and are 
related through the Airy stress function^ y, such that 

0 0 n 

Qy2^ ^yy 92.2 > ^^y g^-Qy- 1 

It follows that the time-independent terms of the stress ten¬ 
sor will vanish in Equation lA.ll and IA.21 but not in Equa¬ 
tion IA.3I The equations of motion for the time-dependent 
displacements thus become 


pui - dx^lx - dyaly = 0 

(A.9) 

pvi - dxaly - dyaly = 0 

(A.IO) 

piv -k kA‘^{wo -k Wi) — 

dx [i<Jxx + <^lx)dx{wo + Wi) 

+ (0'xy + ^ly)9y{wQ + Wl)] - 
9y [{<^xy + ^ly)9x{wo + Wi) 

+ i(^yy + 0-yy)5y(wo -k Wi)] = 0. 

(A.ll) 


When solving these equations, any terms that are not lin¬ 
ear in the derivatives of ui{x,t), vi(x,t) or wi(x,t) can be 
ignored due to small vibrational amplitudes. 

Finite-difference time-domain methods have been used to 
solve Equations IA.HIA.3| As in our previous paper—, the 
equations have been discretized using standard discretiza¬ 
tion schem es^ wit h step sizes Ax = Ay = 0.05 nm and 
At = 0.4y^ dx^/AK = 0.8 fs. The Lame parameters, bend¬ 
ing rigidity and density have been set to the values given 
by the modified Tersoff potential, i.e., y = 167 N m“^, 
A = 23 N m-i, K = 2.8 X lO’^® J and p = 7.42 x 10"'^ 
kg m“®. Fixed boundary conditions are applied in the x 
direction and periodic boundary conditions are applied in 
the y direction. The initial conditions are 


wi(x, 2 /, t = 0) = Re 


(A.12) 


L k 


dtwi {x, y,t = 0) = Re 




with 


Ok = (A.13) 


As in the MD simulations, k. = k^x + kyij is a wavevector 
allowed by the boundary conditions, R = xx + yy is a 
position, A = 0.01 nm is an amplitude and t] = 2 nm is 
the width of the wavepacket. The wavepacket is centered 
around Ro in real space and ko = kxox + kyjj in reciprocal 
space, and a;(fco) is the frequency of out-of-plane vibrations 
with wavevector fcg. 

The static out-of-plane displacement is set to 

u;o(a;, 2 /) = ^1-k asin (A.14) 

where Ly is the system size in the y direction. Fitting to 
the shape of the buckling of the 9.4° boundary produced 
by MD simulations gives Ab = 0.55 nm, ^ = 0.72 nm and 
a = 0.01. As in the MD simulations Ly = 4.5 nm, so m 
must be set to 3 to obtain the correct periodicity in y. The 
system length in the x direction, is set to 100 nm. 

In addition to the static out-of-plane displacement, the 
time-independent terms in the stress tensor components 
are also needed. These have been obtained by fitting to 
the (approximate) stress tensor components obtained from 
MD. Starting with it is seen that if we set 


al, = 


Sin 



(A.I5) 


we obtain a good qualitative correspondence to the MD 
data (see Figure |9]). 

To satisfy the relations between the stress tensor compo¬ 
nents given by Equation lA. 81 we must then set 


xy 


Ay 

27rm 


2 


4 



sin 



^ -2x'^te 

e 


cos 



/ 27Tmy\ 

\^) 

(A.I6) 

(A.17) 


As can be seen in Figures (TU] and [TTJ this functional form of 
the stress tensor components does reproduce the MD result 
in a qualitiative manner. 

In order to compare the results of the continuum mechan¬ 
ical model to those obtained from MD, the transmission 
coefficient Tc is calculated as 


T, = 


/ ^xpY,^.^Q{u)lul{x^,yjAn)+^lvl{x^,yjAn)+^^lwl{x^,yJ,tn)) 


\ 


2^tot 


(A.I8) 


where Wa,, ujy and are the frequencies of vibrations in 
the x, y and z directions, Xi = iAx and yj = jAy indi¬ 
cate a point on the discretization grid, and = nAt is 
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the timestep. The time average is taken over times after 
scattering against the static out-of-plane displacement and 
the total energy is given by 





























Figure 9. (Color online) The stress tensor component cr^j, close 
to the grain boundary (a) as obtained from MD and (b) as 
approximated according to Equation lA. 151 Note that the grain 
boundary is located at a: = 129 nm in the MD simulatons and 
at X = 0 in the continuum mechanical model. 


^tot 


lS.xp 


Lxj'^ 

'y ^ + ^y'^1 Vo ) " 1 “ 

Xi = -L^/2 


(A.19) 


ACKNOWLEDGEMENTS 

The authors would like to thank Prof. Jari Kinaret for 
rewarding discussions. We also acknowledge financial sup¬ 


port from the Swedish Research Council (VR) and the EU 
Graphene Flagship (grant no. 604391). 


^ O. V. Yazyev and Y. P. Chen, Nature Nanotechnology 9 , 755 ^ O. V. Yazyev, Solid State Commun 152 , 1431 (2012). 

(2014). 


7 

























Figure 10. (Color online) The stress tensor component ayy (a) 
as obtained from MD and (b) as approximated according to 
Equation IA.16I 


Figure 11. (Color online) The stress tensor component a^y (a) 
as obtained from MD and (b) as approximated according to 
Equation IA.17I 
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